MODULE OCEAN_ML_DRIVER_V2_MOD

CONTAINS

SUBROUTINE OCEAN_ML_DRIVER_V2 &
 & ( KIDIA    ,KFDIA    ,KLON     ,KLEVO    ,PTSPHY   ,&
 &   PGEMU    ,PSSRFL   ,PSLRFL   ,PAHFS    ,&
 &   PEVAP    ,PUSTOKES ,PVSTOKES ,&
 &   PTAUOCX  ,PTAUOCY  ,PPHIOC   ,PWSEMEAN ,PWSFMEAN ,&
 &   PUO0     ,PVO0     ,PTO0     ,PSO0     ,&
 &   YDCST    ,YDEXC    ,YDMLM    ,&
 &   POTKE    ,PUOE1    ,PVOE1    ,PTOE1    ,PSOE1    )
!
! (C) Copyright 2011- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
! -------
!   THIS ROUTINE IS THE MAIN DRIVER FOR THE OCEAN MIXED LAYER MODEL.

!     ORIGINAL P. A. E. M. JANSSEN   E.C.M.W.F.         21-07-2011
!     BASED ON INFRASTRUCTURE DEVELOPED BY Y. TAKAYA
!
!     PURPOSE
!     -------

!     SOLVES FOR OCEAN SURFACE TEMPERATURE IN THE SURFACE LAYER


!     INTERFACE
!     ---------

!     *OCEAN_ML_DRIVER_V2 IS CALLED BY *SURFPP*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLON*         NUMBER OF GRID POINTS PER PACKET
!     *KLEVO*        NUMBER OF VERTICAL LAYERS

!     INPUT PARAMETERS (REAL):

!     *PTSPHY*       PHYSICS TIME STEP
!     *PGEMU*        THE SINE OF LATITUDE
!     *PSSRFL*       NET SOLAR RADIATION AT THE SURFACE
!     *PSLRFL*       NET THERMAL RADIATION AT THE SURFACE
!     *PAHFS*        SURFACE SENSIBLE HEAT FLUX
!     *PEVAP*        MOISTURE FLUX
!     *PUSTOKES*     SURFACE STOKES VELOCITY X-DIRECTION
!     *PVSTOKES*     SURFACE STOKES VELOCITY Y-DIRECTION
!     *PTAUOCX*      SURFACE MOMENTUM FLUX TO OCEANS X-DIRECTION
!     *PTAUOCY*      SURFACE MOMENTUM FLUX TO OCEANS Y-DIRECTION
!     *PPHIOC*       ENERGY FLUX FROM WAVES TO OCEANS (DIMENSIONAL)
!     *PWSEMEAN*     WINDSEA VARIANCE
!     *PWSFMEAN*     WINDSEA MEAN FREQUENCY


!     *PUO0*         X-COMPONENT OF OCEAN VELOCITY
!     *PVO0*         Y-COMPONENT OF OCEAN VELOCITY
!     *PTO0*         TEMPERATURE PROFILE !!! in degree C !!!
!     *PSO0*         SALINITY PROFILE

!     OUTPUT PARAMETERS (REAL):
!     *POTKE*        UPDATED TURBULENT KINETIC ENERGY IN THE OCEAN MIXED LAYER
!     *PUOE1*        TENDENCY X-COMPONENT OF OCEAN VELOCITY IN TIMESTEP PTSPHY
!     *PVOE1*        TENDENCY Y-COMPONENT OF OCEAN VELOCITY
!     *PTOE1*        TENDENCY TEMPERATURE PROFILE
!     *PSOE1*        TENDENCY SALINITY PROFILE (not used)

 
!     METHOD
!     ------
!     THE MIXED LAYER MODEL HAS AS PROGNOSTIC VARIABLES THE TEMPERATURE 
!     AND CURRENT AT KLEVO LEVELS. THE EQUATIONS FOR MOMENTUM AND 
!     TEMPERATURE ARE SOLVED WITH A SEMI-IMPLICIT SCHEME. SEA STATE 
!     EFFECTS, SUCH AS MIXING BY WAVE BREAKING AND LANGMUIR TURBULENCE AND ALSO 
!     STOKES-CORIOLIS FORCING ARE INCLUDED FOLLOWING THE WORK OF 
!     JANSSEN (2010) BY OBTAINING THE TURBULENT VELOCITY FROM THE SOLUTION OF 
!     THE KINETIC ENERGY BALANCE EQUATION.
! 

! EXTERNALS 
! ---------
!     OC_MLM
!     

! REFERENCE
! ---------
!   P A E M JANSSEN, OCEAN WAVE EFFECTS ON THE DAILY CYCLE IN SST,
!   ECMWF TECHNICAL MEMORANDUM, 2010

! MODIFICATIONS
! -------------
!   
! 
!
!---------------------------------------------------------------------

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST  , ONLY : TCST
USE YOS_EXC  , ONLY : TEXC
USE YOS_MLM  , ONLY : TMLM
USE OC_MLM_MOD   

IMPLICIT NONE

INTEGER(KIND = JPIM), INTENT(IN) :: KIDIA
INTEGER(KIND = JPIM), INTENT(IN) :: KFDIA
INTEGER(KIND = JPIM), INTENT(IN) :: KLON
INTEGER(KIND = JPIM), INTENT(IN) :: KLEVO

REAL(KIND = JPRB), INTENT(IN)    :: PTSPHY 
REAL(KIND = JPRB), INTENT(IN)    :: PGEMU(KLON)  

REAL(KIND = JPRB), INTENT(IN)    :: PSSRFL(KLON) 
REAL(KIND = JPRB), INTENT(IN)    :: PSLRFL(KLON) 
REAL(KIND = JPRB), INTENT(IN)    :: PAHFS(KLON) 
REAL(KIND = JPRB), INTENT(IN)    :: PEVAP(KLON) 
REAL(KIND = JPRB), INTENT(IN)    :: PUSTOKES(KLON)
REAL(KIND = JPRB), INTENT(IN)    :: PVSTOKES(KLON)
REAL(KIND = JPRB), INTENT(IN)    :: PTAUOCX(KLON)
REAL(KIND = JPRB), INTENT(IN)    :: PTAUOCY(KLON)
REAL(KIND = JPRB), INTENT(IN)    :: PPHIOC(KLON)
REAL(KIND = JPRB), INTENT(IN)    :: PWSEMEAN(KLON)
REAL(KIND = JPRB), INTENT(IN)    :: PWSFMEAN(KLON)

REAL(KIND = JPRB), INTENT(IN)    :: PUO0(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(IN)    :: PVO0(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(IN)    :: PTO0(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(IN)    :: PSO0(KLON,KLEVO+1) 

TYPE(TCST),        INTENT(IN)    :: YDCST
TYPE(TEXC),        INTENT(IN)    :: YDEXC
TYPE(TMLM),        INTENT(IN)    :: YDMLM

REAL(KIND = JPRB), INTENT(INOUT) :: POTKE(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(OUT)   :: PUOE1(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(OUT)   :: PVOE1(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(OUT)   :: PTOE1(KLON,KLEVO+1) 
REAL(KIND = JPRB), INTENT(OUT)   :: PSOE1(KLON,KLEVO+1)


INTEGER(KIND=JPIM), PARAMETER :: NSTEP = 3          ! number of sub-iterations
INTEGER(KIND=JPIM) :: JL,JZ

REAL(KIND=JPRB), PARAMETER :: ZEPUST = 0.0001_JPRB  !    security constant for velocity      (m/s)
REAL(KIND=JPRB), PARAMETER :: ZEPS = 1.E-5_JPRB     !    security constant for wave variance
REAL(KIND=JPRB), PARAMETER :: ALPHAD = 100.0_JPRB   !    Defaults nondimensional energy flux into the ocean (water side)

REAL(KIND=JPRB) :: ZDELT, ZWST2, ZROMEGA2, ZZPI, ZINVPTSPHY
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

REAL(KIND=JPRB), DIMENSION(KLON) :: ZPHIOC
REAL(KIND=JPRB), DIMENSION(KLON) :: ZQ_0, ZPHIW, ZWST, ZF, ZZ0
REAL(KIND=JPRB), DIMENSION(KLON) :: ZHS, ZKS, ZU_BC, ZV_BC, PT_BC
REAL(KIND=JPRB), DIMENSION(KLON,KLEVO+1) :: ZUO0, ZVO0, ZTO0, ZSO0

!--------------------------------------------------------------------

IF (LHOOK) CALL DR_HOOK('OCEAN_ML_DRIVER_MOD:OCEAN_ML_DRIVER_V2',0,ZHOOK_HANDLE)
ASSOCIATE(RG=>YDCST%RG, RLVTT=>YDCST%RLVTT, ROMEGA=>YDCST%ROMEGA, &
 & RPI=>YDCST%RPI, &
 & LEOCCO=>YDEXC%LEOCCO, LEOCLA=>YDEXC%LEOCLA, LEOCWA=>YDEXC%LEOCWA, &
 & RHO=>YDMLM%RHO, &
 & RKAP=>YDEXC%RKAP)

! 1. INITIAL SETTING
!    ---------------

ZZPI = 2.0_JPRB*RPI
ZROMEGA2 = 2.0_JPRB*ROMEGA

ZDELT = PTSPHY/REAL(NSTEP)
ZINVPTSPHY = 1._JPRB/PTSPHY


!     2. GENERATE INPUT FIELDS 
   
  DO JL=KIDIA,KFDIA

!  The temperature at the bottom of the profile is imposed to be the input SST (see updclie)
!  !!! except that it is now in degC
!  Impose it as the boundary condition
   PT_BC(JL) = PTO0(JL,KLEVO+1)

!   WATER SIDE FRICTION VELOCITY AND STRESS DIRECTION
    ZWST2    = SQRT(PTAUOCX(JL)**2+PTAUOCY(JL)**2)/RHO
    ZWST(JL) = MAX(SQRT(ZWST2),ZEPUST)

    ZPHIW(JL) = ATAN2(PTAUOCX(JL),PTAUOCY(JL))

!     ENERGY FLUX (defined as negative)
    IF ( PPHIOC(JL) < 0.0_JPRB ) THEN
      ZPHIOC(JL) = PPHIOC(JL)
    ELSE
      ZPHIOC(JL) = -ALPHAD * ZWST(JL)**3
    ENDIF

!     TOTAL HEATFLUX

    ZQ_0(JL) = MAX(PSLRFL(JL)+PAHFS(JL)+RLVTT*PEVAP(JL),ZEPUST)

!     CORIOLIS PARAMETER
 
    ZF(JL) = ZROMEGA2*PGEMU(JL)

!    SIGNIFICANT WAVE HEIGHT AND GRADIENT SCALE STOKES DRIFT

    ZHS(JL) = 4._JPRB*SQRT(MAX(PWSEMEAN(JL),ZEPS))
    ZKS(JL) = MAX((ZZPI*PWSFMEAN(JL))**2/RG,ZEPUST)

!    GRADIENT SCALE PRESSURE TRANSPORT TERM ('ROUGHNESS' LENGTH)

    ZZ0(JL) = 0.5_JPRB*ZHS(JL)
  ENDDO

!
!     3. TEMPERATURE AND CURRENT INCREMENTS DURING ONE PHYSICS 
!        TIME STEP DUE TO SOLAR HEATING, HEATFLUX FROM/TO THE ATMOSPHERE
!        AND UPPER OCEAN MIXING (JANSSEN, 2010)
!
 DO JZ=1,KLEVO+1
   DO JL = KIDIA,KFDIA
     ZUO0(JL,JZ) = PUO0(JL,JZ)
     ZVO0(JL,JZ) = PVO0(JL,JZ)
     ZTO0(JL,JZ) = PTO0(JL,JZ)
   ENDDO
 ENDDO

! if we get currents from ocean model, we could connect them here !!!
 DO JL = KIDIA,KFDIA
   ZU_BC(JL) = 0._JPRB
   ZV_BC(JL) = 0._JPRB
 ENDDO
   
  CALL OC_MLM(KIDIA,KFDIA,KLON,KLEVO,ZDELT,NSTEP,PSSRFL,ZQ_0,ZF,&
 &            ZPHIW,ZHS,ZZ0,ZPHIOC,ZWST,&
 &            PUSTOKES ,PVSTOKES ,&
 &            ZKS,ZU_BC,ZV_BC,PT_BC,&
 &            POTKE,ZUO0,ZVO0,ZTO0,YDMLM,YDCST)
 

 DO JZ=1,KLEVO+1
   DO JL = KIDIA,KFDIA

     PUOE1(JL,JZ) = ZINVPTSPHY*(ZUO0(JL,JZ) - PUO0(JL,JZ))
     PVOE1(JL,JZ) = ZINVPTSPHY*(ZVO0(JL,JZ) - PVO0(JL,JZ))
     PTOE1(JL,JZ) = ZINVPTSPHY*(ZTO0(JL,JZ) - PTO0(JL,JZ))
     PSOE1(JL,JZ) = 0.0_JPRB

   ENDDO
 ENDDO 

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('OCEAN_ML_DRIVER_MOD:OCEAN_ML_DRIVER',1,ZHOOK_HANDLE)

END SUBROUTINE OCEAN_ML_DRIVER_V2

END MODULE OCEAN_ML_DRIVER_V2_MOD
